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Abstract 


A  framework  for  describing  the  deformation  and  failure  responses  of  multi-phase  polycrystalline  microstructures  is 
developed  from  micromechanical  considerations  and  volume  averaging  techniques.  Contributions  from  damage  (i.e.,  dis¬ 
placement  discontinuities  such  as  cracks,  voids,  and  shear  bands)  are  captured  explicitly  in  the  framework’s  kinematics  and 
balance  relations  through  additive  decompositions  of  the  total  deformation  gradient  and  nominal  stress,  respectively. 
These  additive  decompositions — which  notably  enable  description  of  arbitrarily  anisotropic  deformations  and  stresses 
induced  by  damage — are  derived  following  the  generalized  theorem  of  Gauss,  i.e.,  a  version  of  the  divergence  theorem 
of  vector  calculus.  A  specific  rendition  of  the  general  framework  is  applied  to  study  the  response  of  a  dual-phase  tungsten 
(W)  alloy  consisting  of  relatively  stiff  pure  W  grains  embedded  in  a  more  ductile  metallic  binder  material.  In  the  present 
implementation,  a  Taylor  scheme  is  invoked  to  average  grain  responses  within  each  phase,  with  the  local  behavior  of  indi¬ 
vidual  grains  modeled  with  finite  deformation  crystal  plasticity  theory.  The  framework  distinguishes  between  the  effects  of 
intergranular  damage  at  grain  and  phase  boundaries  and  transgranular  damage  (e.g.,  cleavage  fracture  of  individual  crys¬ 
tals).  Strength  reduction  is  induced  by  the  evolving  volume  fraction  of  damage  (i.e.,  porosity)  and  microcrack  densities. 
Model  predictions  are  compared  with  experimental  data  and  observations  for  the  W  alloy  subjected  to  various  loading 
conditions. 

Published  by  Elsevier  Ltd. 


1.  Introduction 

Constitutive  descriptions  for  deterioration  of  material  strength  capacity  due  to  separation  or  rupture  of 
material  have  been  the  focus  of  numerous  investigations  within  the  context  of  continuum  damage  mechanics 
[1].  For  ductile  poly  crystalline  metals,  scalar  damage  descriptions  measuring  porosity  and  reflecting  inelastic 
volume  changes  have  received  a  great  deal  of  attention  in  the  literature  [2,3].  Concurrently,  representations 
based  on  effective  configurations  with  reduced  material  strength  have  been  popularized,  including  models  fea¬ 
turing  scalar  damage  variables  [4]  or  vector-  or  higher-rank  tensor-based  damage  variables  [5,6].  Tensor-based 
treatments  have  also  been  applied  to  describe  degraded  composite  materials  exhibiting  a  nominally  elastic  or 
viscoelastic  response  [7].  In  brittle  ceramics,  scalar  damage  variables  are  frequently  implemented  [8],  although 
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more  complex  methods  have  been  forwarded  to  account  for  anisotropic  strain  rate  accommodation  due  to  dis¬ 
tributed  microcracking  [9].  Non-local  or  gradient-based  measures  have  also  been  used  for  modeling  damage  in 
ductile  [10]  and  brittle  [11]  systems.  From  a  multiscale  modeling  perspective,  it  was  suggested  in  [12]  that  a 
polycrystal  plasticity  theory  can  capture  the  evolution  of  crystallographic  texture  in  conjunction  with  a  scalar 
porosity  description.  Proposed  in  [13]  is  a  micromechanics-based  model  accounting  for  anisotropic  inelastic 
deformation  due  to  grain  boundary  sliding  and  migration. 

Presently,  empirical  scalar-based  damage  descriptions  remain  the  norm  in  practical  numerical  simulations 
of  impact  and  failure  [8,14].  For  ductile  metals,  the  model  in  [4]  features  a  damage  parameter — the  cumulative 
scalar  plastic  strain  at  failure — whose  instantaneous  value  may  depend  upon  the  strain  rate,  temperature,  devi- 
atoric  stress,  and/or  hydrostatic  pressure.  Macroscopic  experiments  (e.g.,  tension  and  torsion  tests)  are  used  to 
determine  failure  strains  and  calibrate  other  material  constants  entering  the  model.  However,  since  parameters 
are  chosen  primarily  upon  consideration  of  macroscopic  (stress-strain)  data,  the  connection  between  these 
empirical  parameters  and  detailed  elements  of  the  micro  structure  is  not  often  evident.  Furthermore,  failure 
properties  are  generally  isotropic  in  the  sense  that  the  evolution  of  a  directionally  invariant  scalar  parameter 
(e.g.,  effective  plastic  strain  at  failure)  dictates  fracture. 

Even  though  micromechanics-based  treatments  characterizing  anisotropic  damage  have  been  forwarded  in 
the  literature  (cf.  [9]),  these  models  have  yet  to  surpass,  in  practical  or  commercial  applications,  the  empirical 
scalar-based  models  that  typically  feature  fewer  material  parameters  and  require  comparatively  less  effort  to 
implement  in  a  numerical  setting.  However,  in  the  near  future,  detailed  multiscale  models  of  damage  with 
explicit  links  between  micro  structural  properties  (e.g.,  grain  size  distributions,  lattice  orientations,  and  grain 
boundary  character)  and  macroscopic  strength  degradation  are  expected  to  come  to  the  fore  in  numerical  sim¬ 
ulations  of  structures  undergoing  failure,  enabling  design  of  materials  for  enhanced  performance  during  fail¬ 
ure  processes  [15],  e.g.,  energy  absorption  in  vehicular  impact  events  or  ballistic  performance  of  armor  and 
projectiles.  Considering  the  steady  improvement  regarding  computational  capabilities  developed  over  the  past 
decade  for  modeling  defects  in  micro  structures,  along  with  experimental  capabilities  for  acquiring  material 
descriptions  and  associated  properties  at  increasingly  fine  length  scales,  multiscale  micromechanics-based 
approaches  akin  to  that  forwarded  here  appear  increasingly  promising. 

Constructed  in  the  present  work  is  a  macroscopic  description  of  damage  evolution  in  multi-phase  polycrys¬ 
tals  via  direct  averaging  of  micromechanical  solutions.  Following  [16],  the  contributions  from  various  damage 
mechanisms  of  arbitrary  geometry  are  explicitly  accounted  for  in  the  deformation  gradient  decomposition, 
leading  to  a  precise  description  of  the  kinematics  of  anisotropic  damage.  While  the  model  is  more  complex 
than  many  of  the  above-mentioned  scalar-based  empirical  approaches,  the  connection  of  material  parameters 
to  the  micro  structure — such  as  lattice  orientation  and  grain  boundary  content — is  more  immediately 
apparent. 

The  focus  of  the  model  presented  here  is  the  thermomechanical  behavior  of  metallic  poly  crystals.  In  crystal 
plasticity  theory  [17],  slip  system  geometries  are  tracked  explicitly  and  as  a  result  account  for  elastic-plastic 
anisotropy.  However,  damage  evolution  on  individual  crystallographic  planes  has  not  been  emphasized  in 
the  literature  for  metallic  crystals  (although  the  model  of  Espinosa  et  al.  [9]  explicitly  accounts  for  fracture 
on  intrinsically  weak  planes  in  ceramic  crystals).  Failure  of  preferred  lattice  planes  can  dominate  the  response 
of  certain  metallic  systems  such  as  body-centered  cubic  (BCC)  tungsten  [18].  Furthermore,  in  multi-phase 
materials,  preferred  orientations  for  strain  accommodation  due  to  damage  depend  upon  the  grain  boundary 
geometry  [19].  Likewise,  in  materials  exhibiting  preferred  void  shapes  and  arrangements,  the  contribution  to 
strain  anisotropy  from  damage  depends  upon  the  stress  state  and  constraints  imposed  by  the  surrounding 
micro  structure. 

In  Section  2,  the  general  framework,  predicated  upon  explicit  volume  averaging  procedures,  is  described. 
Following  suggestions  in  [20,21],  the  primary  mechanical  variables  upon  which  the  framework  is  built  are 
the  net  deformation  gradient  F  and  net  nominal  stress  S.  In  Section  3,  a  particular  version  of  the  general 
framework  is  developed  to  study  the  response  of  a  two-phase  tungsten  alloy  of  high  interest  for  use  in  defense 
applications.  The  assumption  in  [22]  will  be  used  to  account  for  the  deformation  gradient  distribution  within 
each  phase,  with  individual  grains  modeled  via  crystal  plasticity  theory.  Cleavage  fractures  are  dictated  by  a 
traction-based  criterion  on  intrinsically  weak  crystallographic  planes  in  each  W  grain  [18],  while  intergranular 
decohesion  is  controlled  by  a  stress-  and  temperature-based  model  developed  following  consideration  of 
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previous  physical  and  numerical  experiments  [19,23-26].  In  Section  4,  numerical  implementation  of  the  frame¬ 
work  is  described.  Model  predictions  are  then  discussed  in  Section  5:  the  model  adequately  characterizes  the 
mechanical  response  of  a  W  alloy  under  tensile  loading  at  low  and  high  strain  rates  and  under  shear  loading  at 
low  strain  rates.  Conclusions  follow  in  Section  6. 

The  following  notation  is  used.  Vector  and  tensor  quantities  are  represented  with  boldface  type,  while  sca¬ 
lars  and  individual  tensor  components  are  written  in  italics.  The  index  notation  is  often  used,  following  the 
Einstein  summation  convention  and  distinguishing  between  co variant  (subscript)  and  contravariant  (super¬ 
script)  components.  Juxtaposition  implies  summation  over  two  repeated  adjacent  indices  (e.g., 
(AB)a&  =  AacBcb).  The  dot  (scalar)  product  of  vectors  is  represented  by  the  symbol  (e.g.,  a  b  =  aagabbb , 
with  components  of  the  metric  tensor).  Angled  brackets  denote  a  dual  (scalar)  product  (e.g.,  for 
second-rank  tensors,  (A,B)  =  tr(AB)  =  AabBba).  The  colon  denotes  contraction  over  repeated  pairs  of  indices 
(e.g.,  A  :  B  =  tr(ATB)  =  AabBab  and  C  :  D  =  C*bcdDcd).  The  symbol  “(g)”  represents  the  tensor  (outer)  product 
(e.g.,  (a  0  b )ab  =  aabb). 

2.  Multiscale  framework 


2.7.  Kinematics 


Let  X  and  x  denote  local  fine  scale  reference  and  spatial  coordinates  within  a  volume  element  of  material. 
The  local  deformation  gradient  f  is  then  given  by 


r  =  ^ 
ex’ 


r  =  — 

J.A  qxA  ■ 


(1) 


The  net  deformation  gradient  F  for  the  element  is  then  determined  from  the  motion  of  the  volume  element’s 
external  boundary: 


F  =  — 
V 


H 


x  0  N  d.S\  FaA=—  /  xaNA  dS, 


'A 


(2) 


where  x  in  (2)  are  fine  scale  spatial  coordinates  along  external  surface  S  with  outward  unit  normal  N.  The  vol¬ 
ume  of  the  material  element  in  the  reference  configuration  is  denoted  by  V,  and  coincident  uniform  coordinate 
systems  are  assumed  in  each  configuration.  When  the  material  is  damage-free  (i.e.,  the  displacement  is  contin¬ 
uous  and  differentiable  throughout  V,  implying  that  V  is  simply  connected),  quantities  entering  Eqs.  (1)  and 
(2)  are  related  by  the  generalized  theorem  of  Gauss  [16,20,21],  i.e., 


iAf/K‘«!-U£*wr-U 


Cdv, 


(3) 


meaning  that  F  is  the  exact  volume  average  of  f.  However,  when  internal  surfaces  exist  within  V,  (3)  does  not 
apply.  Instead,  the  volume-averaged  deformation  gradient  F  within  intact  material  satisfies 


F  =  —  /  fdV  =  —  [  x  0  N dS  —  2  [ 

V  Jv  V  Js  V  JA 


Nd S--  I  x < 

V  J A 


)  NdcL4, 


(4) 


where  A  denotes  the  union  of  referential  surface  areas  across  which  the  material  displacement  becomes  discon¬ 
tinuous,  with  corresponding  reference  normal  Nd,  by  convention  directed  outward  from  the  reference  surface 
into  the  material.  Quantities  introduced  in  Eq.  (4)  are  illustrated  in  Fig.  1  for  the  particular  case  of  a  polycrys¬ 
tal  undergoing  intergranular  separation,  with  area  A  consisting  of  two  different  grain  boundary  facets  with 
reference  normal  vectors  and  N d2).  Eqs.  (2)  and  (4)  may  be  combined  as 


f  =  f 


,NdS  =  I 


l,dr3L 


x  0  NdcL4  =  F  +  Fd, 


(5) 


emphasizing  the  additive  contributions  of  bulk  material  (F)  and  internal  surface  discontinuities  (Fd)  to  the  net 
deformation  gradient  F  supported  by  the  volume  element. 
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Fig.  1.  Finite  deformation  of  polycrystal  exhibiting  intergranular  fracture. 


When  the  geometry  of  the  damage  entities  may  be  adequately  described  via  a  fewer  number  of  coordinates 
relative  to  the  geometry  of  the  body,  e.g.,  a  planar  crack  in  a  three-dimensional  body  or  a  shear  line  discon¬ 
tinuity  in  a  two-dimensional  body,  Eq.  (4)  acquires  the  reduced  form  [16,19] 

F  =  ^  /  fdE  = /  x(g)Nd^--^  /  [x]  (g)Ndd^,  (6) 

k  Jv  *  JS  "  JA 

s - V - '  S - V - ' 

F  Fd 

where  the  displacement  jump  [x]  =  x+  -  x“  and  the  normal  covector  Nd  =  Nd+  =  -Nd_  on  opposing  (positive 
and  negative)  faces  of  the  discontinuity.  Please  note  that  Eqs.  (l)-(6)  are  valid  for  any  solid  body,  and  make  no 
specific  assumptions  regarding  the  composition  of  the  material  (e.g.,  crystal  structure). 

2.2.  Stresses  and  balance  relations 

Let  s  denote  the  local  nominal  stress  (the  transpose  of  the  first  Piola-Kirchhoff  stress),  related  to  the  fine 
scale  Cauchy  stress  a  by  sAa  =  jf~bu<Jba,  where  in  Cartesian  coordinates  the  Jacobian  determinant  j  =  det(f). 
The  net  nominal  stress  S  for  the  volume  element  is  defined  by 

S  =  -f  f  X  0  tcLS”,  SAa  =J  [  XAfdS,  (7) 

F  Js  V  Js 

where  t  is  the  traction  vector  per  unit  reference  area  acting  on  external  surface  S,  satisfying  ta  =  SAaNA.  The 
conventional  local  balances  of  linear  and  angular  momentum  are  written  as 

—  +  ba  =  Poi?,  (8) 

where  ba  are  components  of  the  body  force  vector  per  unit  reference  volume,  and  po  is  the  reference  mass  den¬ 
sity.  When  the  material  element  is  damage  free  (simply  connected  with  no  internal  surface  discontinuities),  and 
for  the  particular  case  of  locally  quasi-static  conditions  and  in  the  absence  of  local  body  forces,  i.e., 
xa  =  ba  =  0,  there  results  [20,21] 

=  f  pVN.tS  =  i  l  d  V='-Jy  +  sy)  dF  =  i  js<-  dF.  (9) 

However,  when  the  volume  element  contains  internal  surfaces  with  total  area  A,  as  well  as  local  inertia  and/or 
body  forces,  the  nominal  stress  supported  by  the  bulk  material  S  becomes  [19] 

SAa=y  J(^XA+Sya)dV  =  y^XAfdS-y^XA(t‘i)adA=y  Jo*  dV  +  y  j ^  ~  b“)XA  dV 

V _ ^ _ /  V _  v _ ✓  V _ ^ _ /  s _ ^ _ / 

SAa  {Sd)Aa  SAa  (Sb)Aa 


(10) 
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where  S  is  the  volume-averaged  nominal  stress,  Sb  is  the  stress  contribution  due  to  micro-inertia  and  body 
forces,  and  Sd  accounts  for  traction  supported  by  internal  surfaces,  such  that  t  =  Nd  •  s  on  local  internal  sur¬ 
face  element  dA.  Combining  Eqs.  (7)  and  (10),  the  additive  nominal  stress  decomposition  is  obtained: 

S  =  S  +  Sb  +  Sd.  (11) 


Analogously  to  Eq.  (6),  for  a  damage  entity  exhibiting  geometry  of  one  spatial  dimension  less  than  that  of  the 
body,  the  contribution  to  the  stress  from  traction  on  internal  surfaces  can  be  written  as 

Sd  =  i  j(x®[t]cL4,  (12) 

where  X  =  X+  =  X-  and  the  traction  jump  [tj  =  t+  -  t-  =  Nd  •  (s+  -  s~)  on  opposing  (i.e.,  positive  and  neg¬ 
ative)  faces  of  the  internal  surface  A,  for  example  along  opposing  sides  of  a  crack  front.  Notice  that  Sd  van¬ 
ishes  when  all  internal  surfaces  are  traction-free.  At  the  macroscopic  scale  of  observation,  the  stress  balance 
laws  are  here  assumed  to  exhibit  the  same  form  as  Eq.  (8),  i.e., 


a  sAa 

dX1 


+  Ba  =  p0xa , 


' a  riAb 


FaAZ 


'b  r*Aa 


(13) 


where  X  and  x  refer  now  to  macroscopic  referential  and  spatial  coordinates,  respectively,  coincident  with  their 
microscopic  counterparts,  and  with  B  =  V~l  fv  bdV  and  p0  =  E_1  fv  p0dV.  Note  that  Eq.  (13)  are  imposed  by 
assumption  and  are  not  direct  volume  averages  of  Eq.  (8). 


2.3.  Thermodynamics 


The  local  (fine  scale)  balance  of  energy  is  written  as  follows: 

e  +  div0q  —  (s,f)  =  r,  (14) 

with  e  the  internal  energy  per  unit  reference  volume,  q  the  heat  flux  vector  per  unit  reference  area,  and  r  the 
energy  source  per  unit  reference  volume.  Here,  divo  denotes  divergence  with  respect  to  reference  coordinates 
X.  The  local  entropy  inequality  is  written  as  follows,  with  rj  the  time  rate  of  entropy  production  per  unit  ref¬ 
erence  volume  and  0  the  local  temperature: 


The  Helmholtz  free  energy  per  unit  reference  volume  xjj  is  introduced  as 

i/i  =  e  —  Or],  (16) 


from  which,  upon  substitution  of  Eqs.  (14)  and  (16)  into  (15),  the  entropy  relation  becomes 


<M> 


(q,Vog) 

0 


>  i jf  +  0rj, 


(17) 


with  V0  the  covariant  derivative  with  respect  to  X.  The  local  heat  flux  is  dictated  by 

q  =  — k  •  Vo0,  (18) 

where  the  contravariant  conductivity  tensor  k  reduces  to  kAB  =  kdAB  for  isotropic  conduction  in  a  Cartesian 
reference  coordinate  system,  with  SAB  Kronecker’s  delta. 

At  the  macroscale,  relations  analogous  to  Eqs.  ( 14)— ( 18)  are  assumed  to  apply,  with  F  and  S  serving  as  the 
primary  mechanical  variables.  Define  the  following  quantities: 

E=V~l  j  e dV,  ~e=V~x  [  OdV,  R  =  V~'  f  rdV.  (19) 

Jv  Jv  Jv 


The  macrolevel  balance  of  energy  is  then  postulated  as 


E  +  Div0Q  -  (S,  F)  =  R, 


(20) 
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where  Div 0  denotes  divergence  with  respect  to  macroscopic  referential  coordinates,  and  where  the  macroscopic 
heat  flux  Q  satisfies 

Q  =  -K  -  Vo0,  (21) 

with  Vo  the  macroscopic  referential  covariant  derivative  and  K  the  effective  thermal  conductivity  (note  that  K 
is  not  necessarily  a  volume  average  of  k  for  heterogeneous  microstructures).  The  macroscopic  volumetric  free 
energy  is  defined  by 

^  =  E-9fj,  (22) 

and  the  reduced  entropy  equality  is  introduced  at  the  macro-level: 

(?)+!■  (23) 

Here,  fj  is  the  macroscopic  effective  entropy  per  reference  volume  for  the  aggregate,  and  is  not  necessarily  a 
volume  average.  Insertion  of  Eqs.  (20)  and  (22)  into  (23)  then  gives 

(S,F)-<Qfo5>  >'F  +  k,  (24) 

a  relationship  analogous  to  fine  scale  inequality  Eq.  (17).  Please  note  that  Eq.  (24)  is  a  strict  volume  average  of 
Eq.  (17)  under  the  following  conditions:  (i)  damage  and  micro -inertial  effects  are  absent  such  that  F  =  F  and 
S  =  S,  (ii)  isothermal  conditions  exist  such  that  temperature  rates  and  temperature  gradients  vanish,  and  (iii) 
particular  boundary  conditions  are  applied  to  the  surface  S  of  element  V  such  that  the  right  side  of  the  fol¬ 
lowing  equation  vanishes  [20]: 

V~l  J  fsdF  —  FS  =  V~l  J (x  -  FX)  0  n(s  -  S)  dS.  (25) 

Expanding  F  and  S  via  Eqs.  (5)  and  (11),  respectively,  Eq.  (24)  becomes 

(S  +  Sb  +  Sd,F)  +  (S  +  Sb  +  Sd,Fd)  -  v°^  ^  'P  +  dfj.  (26) 

6 


3.  Multiscale  theory  for  dual-phase  poly  crystals 

Here  the  model  framework  of  Section  2  is  applied  to  study  cumulative  deformation  and  damage  in  tungsten 
alloys.  The  particular  material  of  interest  consists  of  relatively  stiff  and  brittle  pure  tungsten  grains  (BCC) 
embedded  in  a  relatively  compliant  and  ductile  matrix  (FCC)  consisting  of  nickel  (50  wt.%),  iron  (25  wt.%), 
and  tungsten  (25  wt.%).  The  composite  microstructure  nominally  is  comprised  of  90%  pure  W  and  10%  matrix 
alloy,  and  thus  features  a  net  weight  distribution  of  93W-5Ni-2Fe.  Typically,  grain  sizes  span  10-30  pm  for 
the  W  crystals  and  200-500  pm  for  the  FCC  phase  [19],  meaning  that  multiple  W  crystals  are  often  embedded 
within  each  single  crystal  of  the  binder  phase.  Because  of  their  relatively  large  mass  density,  high  melting 
point,  and  high  strength  at  elevated  rates  of  loading  [23-26],  tungsten  heavy  alloys  are  attractive  materials 
for  use  in  kinetic  energy  penetrators  (i.e.,  projectiles). 

Under  high  rate  impact  conditions,  tungsten  alloys  may  exhibit  a  complex  set  of  deformation  and  damage 
modes,  each  describable  by  the  model  framework  postulated  in  the  present  paper.  Phenomena  of  primary 
interest  include: 

•  Effects  of  crystallographic  texture  and  grain  elongation.  In  kinetic  energy  penetrators,  swaging  and  pre¬ 
twisting  of  tungsten  rods  may  produce  a  preferred  orientation  of  individual  W  crystals  comprising  the 
micro  structure.  Furthermore,  experimental  and  numerical  studies  have  demonstrated  a  possible  correlation 
between  performance  of  penetrators  and  the  orientations  of  crystals  comprising  their  microstructures 
[27]. 
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Fig.  2.  Finite  tensile  deformation  of  dual-phase  polycrystalline  alloy  exhibiting  transgranular  and  intergranular  fracture. 


Fig.  3.  Intergranular  damage  at  W-W  interfaces  in  93W-5Ni-2Fe  alloy. 


•  Effects  of  a  variety  of  fracture  mechanisms,  as  illustrated  in  Fig.  2.  Under  positive  hydrostatic  stress,  W 
alloys  are  prone  to  failure  by  one  or  more  of  the  following  modes  [19,25,28]:  (i)  intergranular  fracture  at 
W-W  boundaries  (Fig.  3),  (ii)  intergranular  fracture  at  W-binder  phase  interfaces,  (iii)  cleavage  fracture 
of  W  grains,  and  (iv)  intragranular  rupture  of  the  binder  (i.e.,  matrix)  phase.  In  penetration  applications, 
tensile  stress  states  often  arise  when  tungsten  rods  encounter  a  oblique  targets  or  ricochet,  leading  to  bend¬ 
ing  and  subsequent  transverse  fracture.  Preference  of  one  damage  mechanism  over  the  others  may  depend 
on  processing  history,  temperature,  and  loading  rate  [29]. 

•  Effects  of  initial  microstructure:  phase  volume  fractions  and  grain  boundary  contiguity.  The  volume 
fraction  of  binder  phase,  as  well  as  the  connectivity  of  grain  and  phase  boundaries  which  act  as  potential 
initiation  sites  for  fracture,  are  known  to  influence  the  macroscopic  strength  and  ductility  of  the  alloy 
[24-26,29]. 

In  what  follows,  constitutive  models  for  each  relevant  aspect  of  tungsten  behavior  are  addressed  within  the 
context  of  the  framework  of  Section  2.  Specifically,  crystal  plasticity  models  for  tungsten  and  binder  phases, 
intragranular  fracture  criteria,  grain  and  phase  interaction  laws,  and  intergranular  failure  models  are  devel¬ 
oped,  coupled,  and  merged  into  our  framework. 

3.1.  Crystal  plasticity  models  for  microscopic  response 

Crystal  plasticity  models  are  invoked  for  the  response  of  each  crystal  in  the  dual-phase  aggregate.  The 
structure  of  the  model  for  the  tungsten  (BCC)  and  matrix  (FCC)  phases  is  similar,  though  different  parameters 
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are  naturally  used  for  each.  Only  the  essential  elements  are  given  here;  for  additional  details  on  model  devel¬ 
opment,  the  reader  is  referred  to  [19]. 

The  local  deformation  gradient  f  of  Eq.  (1)  is  decomposed  multiplicatively  as 

f  =  ffV,  (27) 

where  P,  f°,  and  f1’  represent,  respectively,  the  kinematics  of  elasticity  and  rigid-body  rotation,  thermal  expan- 
sion  or  contraction,  and  the  cumulative  contribution  of  moving  crystal  defects  (i.e.,  dislocation  glide  and  pseu¬ 
do-slip  due  to  twinning).  The  elastic  and  thermal  terms  dictate  the  deformation  of  the  slip  direction 
contravariant  vectors  s(a)  and  slip  plane  normal  covariant  vectors  m(a): 

s(a)  =  m(«)  =  mfp-'p-1.  (28) 

The  spatial  velocity  gradient  1  is  decomposed  as 
dx  ... 

]  _  _  jre^e-1  fejr0^0-l jre-1  ,  je^pp-l^-l^e-1  (29) 

~  0X  v - - - '  " - v - J 

=le  =\0  =P 

with  the  superposed  dot  the  material  time  derivative.  Of  concern  here  are  cubic  lattices.  Hence,  the  thermal 
deformation  is  assumed  isotropic: 

f  =  fOf-1  =  otT01,  (30) 

where  aT  is  the  thermal  expansion  coefficient  giving  the  change  in  length  per  unit  current  length  per  unit  incre¬ 
ment  in  6 ,  and  1  is  the  unit  tensor.  The  plastic  velocity  gradient  in  the  intermediate  configuration  (denoted  by  b 
in  Fig.  4)  is  defined  as  in  crystal  plasticity  theory  [17]: 

JP  =  fPfP-i  =  sp  y(a)s<a)  ®  ,  (31) 

a=l 

with  y ^  the  plastic  shearing  rate  on  slip  system  a,  spanning  n  total  active  slip  systems. 

We  assume  a  Helmholtz  free  energy  potential  per  unit  mass,  $,  of  the  form 

po->  =  1?  =  19(ee,0,a,  (32) 

where  the  intermediate  configuration  elastic  strain  2(ee )ajS  =  fe“gabfp  ~  ga/h  with  g aj§  a  metric  tensor  on  b.  The 
symbol  £  denotes  a  dimensionless  scalar  internal  variable  representing  stored  micro-elastic  energy  associated 


Fig.  4.  Deformation  maps  and  configurations  of  local  single  crystalline  volume  element. 
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with  crystal  defects  that  may  impede  shearing  on  each  slip  system  (i.e.,  dislocations).  Let  p  denote  the  mass 
density  in  configuration  b  and  let  the  elastic  second  Piola-Kirchhoff  stress  be  written  (se)^  = 


ffl  lccdabfeb  1/5  =  p^jr,  with  f  =  p/p  and  p  the  local  mass  density  in  configuration  b.  The  resolved  Cauchy 


ap 


stress  on  system  a  is  found  by  T(a)  =  a  :  (gs(a)  < 
in  the  spatial  frame  as  follows  [19]: 


m(oc)).  The  localized  energy  balance  in  Eq.  (14)  can  be  written 


pc9  =  t —  p((0^)  —  +  p9doeetf  :  ee  +  div(k  •  Wx9)  +  —  r 

energy  of  lattice  defects 


temperature 

change 


plastic 

dissipation 


thermoelastic  coupling 


heat 

conduction 


Po 


(33) 


heat 

supply 


In  (33),  the  specific  heat  capacity  c  is  defined  as  p0c  =  dee.  Also,  Vx  denotes  covariant  differentiation  with 
respect  to  fine  scale  spatial  coordinates  x,  and  div  likewise  denotes  divergence  with  respect  to  x. 

When  elastic  strains  are  small,  the  instantaneous  response  of  the  material  is  adequately  described  by  linear 
hyperelasticity  theory.  A  free  energy  potential  per  unit  intermediate  configuration  volume  is  specified  in  this 


case  as 


p#=l-ee:C:e:*  +  l-Knm2+y(0), 


(34) 


where  C  and  p  are  the  fourth  rank  elastic  modulus  tensor  in  configuration  b  and  an  effective  shear  modulus, 
respectively,  and  k  is  a  dimensionless  parameter  that  we  assume  is  independent  of  strain  rate  and  temperature. 
The  function  y(9)  =  —  c9\n(9/9o)  accounts  for  the  purely  thermal  energy,  with  0O  a  reference  temperature  at 
which  y  =  0.  From  partial  differentiation  of  (34),  we  see  that  the  stress  satisfies  the  linear-hyperelastic  relation¬ 
ship  se  =  C  :  ee.  For  an  isotropic  response,  the  elastic  modulus  tensor  is  given  by 

=  m~g«P~gxS  +  n(0)(g°xg*S  +  ga5gh),  (35) 

with  Lame’s  constant  l  and  contravariant  components  of  the  metric  tensor  on  configuration  b.  Please  note 
that  (35)  is  appropriate  as  single  crystalline  tungsten  exhibits  the  unique  property  of  isotropy  with  respect  to 
elastic  constants.  Temperature  dependencies  of  moduli  of  (35)  are  modeled  as  follows  to  first  order  for  pure  W: 


(36) 


deA  =  C6U  +  Cen{d  -  273),  dep  =  C%  +  C*,(0  -  273), 
where  CeU:  C°2)1  Celfl,  and  C°2jl  are  material  parameters. 

A  power-law  viscoplastic  flow  rule  is  invoked  to  model  the  time  rate  of  plastic  deformation  within  each 
crystalline  phase  of  the  dual-phase  material: 


?(«) 


t(a)  =  ro(^)  sgn(xW). 


(37) 


In  (37),  i’o  and  m  are  material  constants,  g>J>  is  the  slip  resistance,  T,a)  =  and  sgn(x)  =  x/\x\,  with 

sgn(O)  =  1.  Thermal  softening  attributed  to  increased  dislocation  mobility  at  high  temperatures  is  incorpo¬ 
rated  via  the  power-law  form 


g 


(38) 

with  C'  the  flow  resistance  at  reference  temperature  do  and  p  a  dimensionless  constant.  The  following  rela- 
tionship  is  imposed  between  the  “average”  hardening  over  all  systems  at  fixed  reference  temperature  and 
the  internal  variable  £: 


\  Llteo0  ~  4a))  =  *Pb\fPri 


(39) 


=4 


with  gW  an  initial  yield  stress,  b  the  magnitude  of  the  Burgers  vector,  and  pT  the  total  dislocation  line  length 
per  unit  intermediate  configuration  volume  associated  with  shearing  impedance.  The  scalar  proportionality 
factor  a  accounts  for  dislocation  interactions.  Both  lattice  friction  stress  and  effects  of  initial  dislocation  den¬ 
sity  are  incorporated  in  the  initial  yield  stress  g 
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For  the  BCC  W  phase,  we  allow  slip  in  the  (111)  close-packed  directions  on  any  of  the  {110}  and  {112} 
families  of  planes,  meaning  the  number  of  potentially  active  slip  systems  is  n  =  24.  Evolution  of  slip  resistance 
at  reference  temperature  60  is  dictated  by  a  hardening-minus-dynamic-recovery  relation: 

£oa)  =  AY  I  “  Bgoa)  Y  If  I ’  (40) 

with  the  interaction  matrix  satisfying 

q;  =  s;+q(  i -«?;),  (4i) 

where  q  is  the  latent  hardening  ratio. 

In  the  W-Ni-Fe  matrix  material,  the  number  of  potentially  active  slip  systems  is  chosen  as  n  =  12.  Assume 
that  dislocations  glide  in  (1 1 0)  close-packed  directions  on  { 1 1 1}  planes  for  this  FCC  metal.  Elastic  isotropy  is 
also  assumed  for  this  phase,  meaning  that  Eq.  (35)  applies,  albeit  with  elastic  stiffness  constants  for  the  matrix 
substantially  lower  in  magnitude  than  those  for  the  pure  W.  We  specify  strain  rate-  and  temperature-depen¬ 
dent  slip  resistances  in  crystals  comprising  the  more  compliant  matrix  phase  via  Eqs.  (40)  and  (41),  though 
with  different  values  of  A ,  B ,  and  q  than  those  invoked  for  the  pure  W  grains. 

Material  properties  for  thermoelasticity  and  plasticity  in  each  phase  are  listed  in  Table  1. 

3.2.  Fracture  model  for  W  cleavage 

Pure  tungsten  single  crystals  are  known  to  fracture  along  the  preferred  crystallographic  planes  having  ori¬ 
entations  {100}  and  {110},  with  the  former  characterized  by  relatively  lower  fracture  toughness  than  the  lat¬ 
ter  [18].  Furthermore,  tungsten  exhibits  a  ductile-to -brittle  transition,  with  cleavage  resistance  increasing 
dramatically  at  ambient  temperatures  above  around  370  K  [18,29].  Despite  the  increase  in  fracture  toughness 
of  W  with  temperature,  O’Donnell  and  Woodward  [29]  recorded  an  increase  with  temperature  in  contributions 
from  W  cleavage  relative  to  the  influence  of  intergranular  fracture  and  matrix  rupture  mechanisms  in  the 
tensile  response  of  the  dual-phase  tungsten  alloy  of  primary  interest  in  the  present  work.  The  fracture  process 
is  thought  to  commence  on  weak  tungsten-tungsten  grain  boundaries,  and  then  proceed  to  subsequent  grain 


Table  1 


Thermoelastic  and  plastic  properties  for  crystalline  phases 


Parameter 

Value  (W) 

Value  (matrix) 

2 

204  Gpa 

137  GPa 

P 

161  Gpa 

99  GPa 

Po 

19350  kg/m3 

9200  kg/m3 

c 

134  J/(kg  K) 

382  J/(kg  K) 

7o 

0.001 

0.001 

m 

20 

20 

q 

1.4 

1.0 

A 

630  MPa 

200  MPa 

B 

1.5 

0.4 

gf 

500  MPa 

150  MPa 

P 

-1.5 

-1.5 

do 

300  K 

300  K 

0CT 

5.3(10)"6/K 

1.5(10)“5/K 

k 

160  W/(m  K) 

100  W/(m  K) 

a 

0.439 

1.03 

b 

0.275  nm 

0.257  nm 

K 

1333 

200 

cu 

-3.4  MPa/K 

- 

ce 

0.0065  MPa/K2 

- 

ce 

-10.3  MPa/K 

- 

ce 

-0.0041  MPa/K2 

- 
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cleavage  when  the  number  of  intergranular  fracture  sites  is  insufficient  to  support  crack  propagation  across  the 
specimen  [19,25,28,29]. 

The  present  model  focuses  upon  tensile  cleavage  fracture  under  mode  I  conditions.  A  stress-  and  temper¬ 
ature-based  initiation  criterion  is  proposed,  along  with  a  micromechanically  inspired  model  for  crack  exten¬ 
sion  on  the  fracture  plane  and  crack  opening  normal  to  the  plane  of  cleavage.  Tangential  initiation  criteria  and 
opening  displacements  are  not  considered  here,  as  the  material  is  known  to  exhibit  much  greater  ductility 
under  shearing  and  compressive  loading  conditions  [23].  This  model  accounts  for  the  microstructural  config¬ 
uration  on  fracture  (through  the  orientation  of  preferential  cleavage  planes),  and  permits  anisotropic  global 
deformation  due  to  damage. 

The  initiation  criterion  for  cleavage  fracture  is  specified  as 


Atf  =  sm(0), 


(42) 


where  sAa  is  the  local  nominal  stress,  mf  and  My  are  spatial  and  referential  co-vectors  normal  to  potential 
cleavage  planes  with  index  k,  and  F\d)  is  a  temperature-dependent  fracture  strength  obeying  the  linear 
relationship 


Ci100}  +  c\loo\0  -  00)  for  {100}  planes, 
C{110}  +  —  do)  for  {110}  planes. 


Values  of  parameters  in  Eq.  (43)  were  selected  upon  assuming  proportionality  between  our  fracture  stress  if® 
and  the  published  fracture  toughness  data  of  pre-cracked  tungsten  single  crystal  bars  loaded  in  three-point 
bending  [18].  Values  are  listed  in  Table  2.  Upon  initiation  of  damage,  i.e.,  satisfaction  of  criterion  in  Eq. 
(42),  crack  growth  and/or  opening  may  occur  on  damage  plane  k.  It  is  assumed  here  that  each  grain  or  sub¬ 
grain  can  support  only  one  microcrack,  i.e.,  fracture  is  limited  to  the  first  plane  achieving  the  criterion  in 
Eq.  (42).  The  contribution  of  cleavage  damage  to  the  global  deformation  gradient  F  is  calculated  from 
Eq.  (6),  where  the  subscript  W  corresponds  to  the  pure  tungsten  phase: 


Fd 

—  — 


Hj 


w = 1 

k= i 


)Ndd^w  =  y^ 

W 


)  N^)W  =  —  V(Ad  (g>  M)W. 


NW 


(44) 


Here  summation  runs  over  7VW  distinct  damage  planes,  each  with  corresponding  referential  area  •  The 
volume  fraction  of  pure  W  in  the  dual-phase  system  satisfies  /w  =  Fw/F,  with  Ew  =  ^  the  total  volume 
of  the  W  crystals.  For  mode  I  cleavage  on  preferred  planes,  the  reference  orientation  satisfies  =  M^, 
and  the  direction  of  opening  is  collinear  with  m^,  i.e.,  =  ||^d^  ||m^.  The  magnitude  of  the  crack  opening 

displacement  jump  (multiplied  here  by  the  microcrack  area)  follows  from  [9],  assuming  circular- shaped 
flaws: 


||Ad«|| 


(  I6n(l  —  v2) 

J  3  E 

0 

^  0 


for  ^  0, 

for  <  0. 


(45) 


Table  2 

W  grain  cleavage  properties 


r{100} 

1.75  GPa 

r{100} 

0.010  GPa/K 

c(110} 

1.90  GPa 

r{110} 

^2 

0.010  GPa/K 

rR 

2.66  km/s 

r° 

bW 

104/s 

zl 

1.0 

z2 

1.0 

jjnax 

rw 

30  pm 
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In  Eq.  (45),  v  and  E  are  Poisson’s  ratio  and  Young’s  modulus  for  undamaged  pure  crystalline  tungsten,  and 
the  crack  opening  displacement  is  non-zero  only  for  tensile  loading,  i.e.,  for  s ^  >  0.  The  microcrack  radius 
and  area  are  related  by  .  The  time  rate  of  damage  deformation  is  thus  dictated  by  the  rate  of 

nominal  stress  as  well  as  the  rate  of  crack  extension,  the  latter  described  by  a  strain  rate-  and  stress-dependent 
growth  law  of  the  form 


r®  - 

'  W  - 


'W 


#) 


Cl(w)  (  1  — )  ^0  for  s ^  >  s ^  and  rfj  < 


,(*) 


w 


rmax 
rW  ’ 


for  ^  s®  or  r$  =  r\ 


w 


(46) 


Above,  is  the  Rayleigh  wave  speed  in  pure  W,  the  effective  strain  rate  s  =  2/3)d  :  d,  with  2 dab  =  lab  + 

lba,  and  is  a  normalization  parameter  required  on  dimensional  grounds.  Eq.  (46)  is  a  stress-  and  strain-rate 
based  generalization  of  the  fracture  toughness-based  approaches  of  [8,9].  However,  in  contrast  to  these  ap¬ 
proaches,  (46)  does  not  require  an  assumption  on  the  initial  crack  size  in  order  to  propagate  damage,  such 
that  one  may  assume  A$t=, 0  =  0  without  difficulty.  The  conditions  for  crack  growth  prevent  extension  when 

the  driving  force  and  ensure  that  the  crack  diameter  does  not  exceed  the  physical  dimension  of  grain 

k ,  i.e.,  the  grain  size  r^ax.  Exponents  Zi  and  z2  control,  respectively,  the  relative  magnitudes  of  strain  rate  and 
stress  influences  on  crack  extension.  Properties  for  the  cleavage  model  applied  to  the  particular  W  heavy  alloy 
considered  here  are  given  in  Table  2. 

The  referential  orientation  unit  vector  in  Eqs.  (42)  and  (44)  remains  invariant  with  time,  and  depends 
only  upon  the  initial  crystallographic  orientation.  In  contrast,  reorientation  of  the  spatial  orientation  vector 
follows  from 


m(*)=reM(*),  (47) 

where  re  is  the  elastic  rotation  of  grain  k  associated  with  the  polar  decomposition  of  the  elastic  deformation 
gradient  f  =  reue.  Thus,  under  a  superposed  global  rigid  body  rotation  such  that  F  — ►  QF,  with  QT  =  Q-1  and 
detQ  =  1,  the  consistent  transformations  are  f  —*  Qf,  re  — *  Qre,  and  F^  —*  QF^. 

A  strain-like  symmetric  tensor  associated  with  damage  in  crystallite  k  is  introduced  as 

(M  \2 

C^4)=  (S)ll[x]llJ  (Nd®N (48) 
the  trace  of  which  is  found  as 


Note  that  for  purely  mode  I  fractures,  where  each  corresponding  displacement  jump  and  normal  vector  to  the 
crack  plane  are  co-linear,  the  volume  fraction  of  damage  per  unit  reference  volume  of  the  local  grain  k  is  then 

found,  in  Cartesian  coordinates,  by  the  quantity 


3.3.  Intergranular  fracture  model 


The  model  discussed  here  collectively  accounts  for  the  following  three  mechanisms  of  damage  in  the  tung¬ 
sten  alloy:  separation  at  W-matrix  interfaces,  failure  along  W-W  grain  boundaries,  and  matrix  rupture  (which 
is  generalized  here  to  include  failure  at  interfaces  between  two  matrix  grains  as  well  as  crack  propagation 
across  matrix  grains).  The  matrix  rupture  mechanism  is  labeled  here  as  an  intergranular  failure  mode,  as 
we  do  not  prescribe  a  damage  criterion  for  cleavage  within  each  individual  matrix  grain. 

Fractures  initiate  on  preferred  locations  in  the  micro  structure.  Experiments  have  implied  that  weakest 
intergranular  failure  sites  in  the  alloy  under  consideration  are  W-W  contacts,  followed  by  W-matrix  phase 
boundaries,  followed  by  internal  boundaries  within  the  matrix  phase  [24,25].  However,  a  previous  numerical 
investigation  [19]  demonstrated  a  tendency  for  initiation  on  W-W  boundaries,  even  if  all  interfaces  are  equally 
strong,  simply  due  to  relatively  larger  stresses  supported  by  the  W  grains.  Relative  strengths  of  the  interfaces 
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are  thought  to  vary  with  temperature  and  processing  conditions  that  may  affect  impurity  concentration  at  the 
interfaces  [28,29]. 

Analogously  to  Eq.  (42),  stress-controlled  damage  initiation  is  assumed.  This  assertion  follows  from  exper¬ 
imental  agreement  in  peak  resolved  normal  failure  stresses  observed  in  macroscopic  tension,  bending,  and 
spallation  [25,28].  Criteria  are  written  in  a  general  form  as 

S(i)  _  (50) 

where  m®  and  My  are  spatial  and  referential  unit  normal  vectors  to  potentially  damaged  planes,  and  the  index 
i  spans  all  potential  fracture  sites  for  the  mechanisms  of  W-W,  W-matrix,  and  matrix-matrix  separation. 
Components  of  the  total  stress  SAa  are  defined  in  Eq.  (7).  Fracture  stresses  S®  are  assumed  linearly  temper¬ 
ature-dependent  and  written  as 

5«(0)  =  CI1+C'(0-0O).  (51) 

Here  C\  and  C\  are  material  parameters  and  9  is  the  average  temperature  introduced  by  definition  (19).  Note 
that  Eq.  (51)  assumes  normal  stress-based  initiation;  shear  stresses  are  not  assumed  to  initiate  damage,  follow¬ 
ing  previous  arguments  and  experimental  observations.  The  number,  initial  orientation  M%\  and  maximum 
size  of  each  potential  damage  plane  must  be  assigned  as  initial  conditions,  based  on  micro  structural  charac¬ 
terization.  In  the  implementation  that  follows,  we  assume  a  distribution  of  damage  planes  exists  encompassing 
all  possible  orientations,  thereby  reducing  Eq.  (50)  to  a  maximum  normal  stress  criterion  with  Ni  =  1.  In  other 
words,  we  assume  the  existence  of  a  single  crack  plane  of  initial  orientation  =  M,  where  M  is  the  principal 
direction  corresponding  to  the  maximum  positive  eigenvalue  of  the  symmetric  nominal  stress  tensor  SAaF~1B. 
Furthermore,  damage  planes  are  assumed  to  rotate  with  the  average  material  deformation  F  of  Eq.  (4)  in 
Cartesian  coordinates  as 

m(f)  =RM(f),  (52) 

where  the  polar  decomposition  of  the  material  deformation  is  written  as  F  =  RU. 

Upon  initiation  of  damage  through  attainment  of  criterion  in  Eq.  (50),  crack  growth  and/or  opening  may 
occur  on  damage  plane  i.  The  contribution  of  damage  to  the  global  deformation  gradient  F  is  calculated  from 
Eq.  (6),  where  the  subscript  label  I  corresponds  to  the  three  aforementioned  intergranular  mechanisms  of 
interest  here: 


Fd 


1 

V 


0Ndd A 


N  i 

£([xH0NO(i) 


1 

V 


N I  _ 

yyid  <g>  M)(i). 

i=i 


(53) 


Note  that  Ai  spans  failure  sites  along  W-W  interfaces,  W-matrix  sites,  and  matrix  internal  boundaries,  and 
that  for  mode  I  fractures,  Nd(z)  =  and  Ad(z)  =  ||Ad(z)||m(/).  From  Eqs.  (3),  (4),  (52)  and  (53),  under  super¬ 
posed  rigid  body  rotation,  F  —>  QF,  f  — ►  Qf,  and  F  — >  QF,  leading  to  the  consistent  relation  Fd  — ►  QFd. 

The  magnitude  of  intergranular  crack  opening  (multiplied  here  by  the  crack  area,  assuming  a  circular  crack 
of  radius  r^),  is  described  by  a  relation  analogous  to  Eq.  (45): 


||Ad«||  =  < 

"  167i(l  —  v2) 

3  E 

5(i)  (r[°)3 

^  0  for  S' !  ^  0, 

(54) 

0 

for  Sil)  <  0. 

Here  v  and  E  are  the  effective  Poisson’s  ratio  and  Young’s  modulus  for  the  undamaged  dual-phase  tungsten 
alloy,  and  the  crack  opening  displacement  is  non-zero  only  for  tensile  loading.  The  microcrack  radius  and  area 

are  related  by  A |z)  =  7i(rjz^  .  The  rate  of  crack  growth  is  again  described  by  a  stress-  and  strain  rate-driven 
growth  law  of  the  form 


Cf  ^  ^1  -  ^  Ss  0  for  St,]  >  S{,)  and  r(  <  r”Yd\ 

k  0  for  S(l>  y  Slr>  and  r(  =  r”ax. 


(55) 
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Table  3 


Intergranular  fracture  properties 


C’l 

1.20  GPa 

c\ 

0.010  GPa/K 

cf 

2.57  km/s 

(2/3)  104/  s 

Z 

1.0 

Z2 

1.0 

E 

366  GPa 

V 

0.29 

^max 

100  jam 

In  Eq.  (56),  Cf  is  the  Rayleigh  wave  speed  in  the  dual-phase  alloy  (computed  from  the  effective  elastic  con¬ 
stants  of  the  undamaged  composite),  the  effective  strain  rate  E  =  (2/3)D  :  D,  with  2 Dab  =  Lab  +  Lba, 

L  =  FF-1,  and  ej,  Z i,  and  Z2  are  material  constants  analogous  to  those  introduced  in  Eq.  (46).  Crack  exten¬ 
sion  is  prevented  when  the  driving  force  S®  ^  and  the  crack  diameter  is  forbidden  from  exceeding  the 
physical  dimension  r^ax  of  oriented  surface  i.  Material  parameters  for  the  intergranular  damage  model  are 
provided  in  Table  3. 

A  strain-like  symmetric  tensor  associated  with  intergranular  damage  in  V,  here  restricted  to  a  single  crack 
plane  (Nj  =  1),  is  introduced  as 

c?=  (yllMIl)  (Nd®Nd),  (56) 

the  trace  of  which  is  found  as  follows,  since  Nd  is  a  unit  vector: 

tr(C?)=(f||[x]||)2.  (57) 

For  purely  mode  I  microcracking,  the  volume  fraction  of  intergranular  damage  per  unit  reference  volume  of 
the  aggregate  is  then  equal  to  yj tr(Cd). 

3.4.  Grain  and  phase  interactions 


The  general  framework  introduced  thus  far  allows  one  to  opt  for  one  of  a  variety  of  interaction  schemes  for 
grains  of  like  and  differing  phases,  for  example  Taylor  [22]  or  self-consistent  schemes  [30],  or  fully  resolved 
calculations  at  the  microscale  in  which  grain  and  phase  topologies  are  fully  described  [19].  The  latter  two 
methods  are  accompanied  by  increased  complexity  with  regards  to  problem  solution  (e.g.,  numerical  imple¬ 
mentation  and  computation  times).  Fully  resolved  micro  structural  calculations,  in  which  individual  integra¬ 
tion  points  encompass  volumes  on  the  scale  of  a  few  micrometers,  are  not  yet  tractable  for  solving  complex 
boundary  value  problems  with  global  dimensions  on  the  order  of  more  than  several  millimeters. 

In  the  present  implementation  we  employ  Taylor’s  assumption  [22]  of  deformation  gradient  uniformity 
among  phases  comprising  each  crystalline  volume  element  V.  The  average  material  deformation  of  Eq.  (4) 
is  rewritten  for  the  particular  alloy  of  study  as 


/w  V^f< 


NW 


F  = 


N' 


i=  1 


r  NM 

J M  X^fU) 


M 


rw 


(58) 


where  subscripts  W  and  M  denote  grains  of  pure  W  phase  and  matrix  binder  phase,  respectively,  N  and  / rep¬ 
resent  the  number  of  grains  of  equal  volume  and  volume  fraction  of  the  particular  phase,  and  f 1 ^  (f  ^ )  is  the 
local  deformation  gradient  assumed  uniform  within  grain  i  (j)  of  the  particular  phase.  The  response  of  each 
grain  is  dictated  by  the  constitutive  theory  of  Section  3.1.  The  deformation  constraint  imposed  at  the  phase 
level  is 
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F=/w'Fw=/m1Fm,  (59) 

and  that  at  the  grain  scale  within  each  phase  is 

/w'Fw  =  ?w>  /M'FM=e  (60) 

Combination  of  Eqs.  (59)  and  (60)  results  in  a  uniform  deformation  gradient  f  for  each  grain  within  V,  regard¬ 
less  of  phase,  leading  to  the  well-known  upper  bound  on  effective  stiffness  for  the  aggregate,  in  the  absence  of 
damage.  Previous  work  [19]  has  indicated  that  under  high  rate  tensile  loading  conditions,  the  ductile  matrix 
phase  tends  to  accommodate  a  larger  magnitude  of  strain  relative  to  the  stiffer  W  phase,  though  at  larger 
deformation  levels,  assumption  Eq.  (59)  may  be  more  justifiable  as  strain  hardening  due  to  dislocation  accu¬ 
mulation  leads  to  a  more  uniform  tangent  stiffness  among  phases  [28].  Regardless,  Eq.  (59)  can  be  used  to 
compare  predictions  of  material  response  for  various  micro  structural  configurations  with  the  caveat  that 
numerical  results  may  be  overly  stiff. 

Average  material  stresses  for  the  aggregate  are  found  as  follows,  specialized  here  to  the  two-phase  tungsten 
alloy  of  interest: 


where  and  are  local  nominal  stresses  in  W  and  matrix  phases  for  each  grain  i  or  j.  Whereas  each  grain  is 
assigned  the  same  deformation  gradient  through  Eq.  (58),  different  stresses  arise  among  grains  due  to  the 
micro  structure  (phase  and  lattice  orientation)  and  corresponding  material  properties.  Note  that  internal 
surfaces  are  presumed  traction  free,  such  that  the  interface-induced  stress  Sd  =  0  in  Eqs.  (10)-(12).  This 
assumption  is  thought  to  be  warranted  under  the  mode  I  fracture  modes  modeled  here,  though  may  not  be 
justified  should  the  model  be  extended  to  account  for  tangential  sliding  and  friction  at  failed  interfaces.  As 
individual  grains  are  not  resolved  spatially,  stresses  due  to  micro-inertia  are  not  computed,  i.e.,  we  set 
Sb  =  0  in  Eq.  (11).  In  the  absence  of  body  forces,  this  assumption  corresponds  to  uniform  acceleration  xa 
of  Eqs.  (8)  and  (10)  throughout  each  grain  in  V  when  referred  to  a  local  referential  coordinate  system  located 
at  the  centroid  of  each  grain.  Macro-inertia  of  Eq.  (13)  is  fully  taken  into  account,  however,  as  the  present 
assumption  corresponds  to  a  uniform  xa  assigned  over  each  V. 

In  the  present  implementation,  thermal  interactions  are  neglected  between  grains  and  phases,  meaning  that, 
depending  upon  the  applied  deformation  rate,  we  impose  in  Eq.  (33)  either  the  isothermal  constraints  6  =  0 
(low  strain  rates  and  long  deformation  histories)  or  adiabatic  constraints  (k,  Vx0)  =  0  (high  strain  rates  and 
short  deformation  histories),  throughout  the  volume  V.  This  assumption  is  necessary  within  the  context  of 
the  Taylor-type  approximations  above,  as  temperature  gradients  between  neighboring  grains  cannot  be  spa¬ 
tially  resolved. 

3.5.  Homogenization  of  damage 

Effects  from  mechanisms  modeled  individually — elastoplasticity  within  each  phase,  cleavage  fracture  of  W 
grains,  and/or  intergranular  damage  modes — are  homogenized  in  a  scheme  that  adheres  to  the  general  frame¬ 
work  of  Section  2.  The  total  deformation  gradient  for  the  element  V  is  computed  from  Eq.  (5)  as 

F  =  F  +  Fd  =  F  +  F^+F?,  (62) 

Fd 

where  F  includes  the  volume-averaged  deformation  gradients  from  tungsten  and  matrix  phases,  and  and  Fd 
are  contributions  from  W  cleavage  and  intergranular  mechanisms  described  by  Eqs.  (44)  and  (53),  respec¬ 
tively.  Total  stresses  entering  Eq.  (11)  are  calculated  as  follows: 

S  =  A  :  S,  SAa  =  AAfSBb ,  (63) 

where  the  rank  four  object  A  accounts  for  stress-  and  tangent  stiffness  reduction  due  to  intergranular  damage. 
Note  that  A  is  needed  in  the  present  implementation  of  the  model,  as  strength  reduction  due  to  damage  is  not 
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accounted  for  in  Eq.  (61)  or  the  Taylor-type  averaging  scheme.  However,  should  the  framework  be  extended 
to  fully  resolved  calculations  at  the  microscale  in  which  damage  entities  at  grain  and  phase  boundaries  are 
described  individually  (e.g.,  [19]),  and  their  effects  on  the  local  stress  fields  within  grains  are  captured  explicitly, 
one  may  set  A  =  1  0  1.  In  the  language  of  continuum  damage  mechanics  [1,5,6],  S  of  Eq.  (63)  may  be  regarded 
as  an  'effective  stress’,  and  A  a  'damage  effect  tensor’. 

As  a  first  approximation,  we  assume  that  the  average  material  nominal  stress  is  reduced  in  an  isotropic 
manner  due  to  microcracking,  i.e., 

<B  =  ASiSl,  (64) 

where  A  is  a  scalar  damage  variable  depending  in  a  multiplicative  fashion  upon  the  volume  fraction  of  damage 
(i.e.,  porosity),  density  of  cleavage  cracks,  and  density  of  intergranular  cracks: 

f  (1  —  OL<p(p)(l  -  (l  -  for  cp>  0, 

^  =  \  1  for  cp  =  (X 

[  0  for  (p  ^  <pfai1, 

The  fraction  of  cumulative  damage  per  unit  reference  volume, 
quantities  introduced  in  Eqs.  (49)  and  (57): 

*=  +  yWfl  • 

v  y  volume  fraction 

volume  fraction  of  transgranular  damage  of  intergranular  damage 


(65) 

coi  ^  oj[ai1  or  cow  ^  co^1- 

cp,  is  defined  by  superposition  of  kinematic 


Scalar  crack  densities  are  defined  in  the  usual  manner  as  effective  cracked  areas  per  unit  reference  volume: 


(67) 


Maximum  crack  densities  supported  by  each  mechanism,  denoted  by  co™x  and  co^ax,  are  found  by  substituting 
respective  micro  structure  parameters  r™x  and  r™ax  from  Eqs.  (46)  and  (55)  into  each  of  Eqs.  (67).  Critical 
porosity  and  crack  densities  above  which  total  failure  (i.e.,  total  stress  relief)  occurs  can  be  specified  in  Eq. 
(65)  by  parameters  <pfai1,  co and  co^11.  Relative  effects  of  porosity,  crack  density  on  W  planes,  and  crack  den¬ 
sity  at  intergranular  sites  are  weighed,  respectively,  by  the  scalar  material  constants  a^,  aw,  and  oq.  Note  that 
in  the  present  implementation,  strength  reduction  from  damage  would  not  occur  under  conditions  of  purely 
compressive  or  shear  loadings,  since  only  mode  I  fractures  are  considered  in  the  fracture  models;  hence,  the 
specification  in  Eq.  (65)  that  A  —  1  when  cp  =  0. 

The  form  of  Eq.  (65)  was  postulated  in  part  upon  consideration  of  the  following  special  cases.  Setting 
a ^  =  1  and  aw  =  oq  =  0  results  in  A  =  1  -  cp,  the  usual  prescription  for  the  reduction  in  effective  stiffness  in 
an  elastic  body  supporting  a  dilute  concentration  of  voids  (cf.  [3]).  Setting  oq  =  1  and  =  aw  =  0  results 
in  A  =  1  —  cdi/cq™3*,  an  appropriate  reduction  factor  for  the  effective  axial  stiffness  (or  stress)  of  a  body  con¬ 
taining  a  single  (intergranular)  flaw  of  radius  oriented  perpendicular  to  the  loading  direction  [5].  In  the  pres¬ 
ent  implementation  wherein  a  mixture  of  damage  mechanisms  evolves,  values  of  parameters  entering  (65), 
listed  in  Table  4,  were  chosen  such  that  the  model  predictions  agree  with  experimental  data  and  observations, 
as  will  be  discussed  more  in  Section  5. 


Table  4 

Effective  damage  properties 


Ucp 

1.00 

aw 

1.00 

OCl 

0.15 

<v<ax 

0.38 

oo 

<pfaii 


oo 
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Thermodynamic  considerations  are  now  revisited  briefly  for  the  particular  model  of  Section  3.  The  macro¬ 
scopic  energy  balance  Eq.  (20)  and  dissipation  inequality  (26)  become,  within  the  context  of  constitutive 
assumptions  Eqs.  (63)  and  (64)  and  adiabatic  conditions, 


(AS,F)  +  (AS,iiw  +  ii1)=±(¥+6rj), 

(68) 

(ASP)  +  (/IS,  F(v  +  F ?)  >  IP  +  ~6rj . 

(69) 

4.  Numerical  implementation 


The  plasticity  and  damage  model  formulated  in  Section  3  was  implemented  numerically  in  order  to  dem¬ 
onstrate  its  predictive  capabilities  and  compare  these  predictions  with  experimental  findings.  In  the  present 
work  the  model  was  implemented  within  a  three-dimensional  material  point  simulator,  which  could  act  as 
a  single  integration  point  within  a  three-dimensional  Lagrangian  finite  element,  for  example.  Deformation- 
controlled  simulations  were  conducted,  with  imposed  values  of  F  (see  Eqs.  (58)-(60))  applied  to  the  grains 
of  the  micro  structure.  Thermal  conditions  were  specified  as  either  isothermal  or  adiabatic. 

For  numerically  integrating  the  thermomechanical  response  history  of  crystals  comprising  pure  W  and  bin¬ 
der  phases  according  to  the  constitutive  model  of  Section  3.1,  implicit  techniques  were  used,  as  discussed  in 
detail  in  [19],  and  described  briefly  in  what  follows.  Let  subscripts  t  and  t  +  At  denote  consecutive  computa¬ 
tion  instances  in  a  nonlinear  analysis,  i.e.,  start  and  end  times  in  a  particular  iteration.  Inelastic  shearing  rates 
for  an  increment  spanning  times  t  and  t  +  At  are  found  implicitly  from  values  of  the  resolved  shear  stress  and 
hardening  variables  at  the  end  of  the  iteration: 


y(a) 


7o 


r(a) 

lt+At 


(a) 

&t+At 


(70) 


An  iterative  procedure  is  invoked  to  solve  Eq.  (70),  as  T^At  and  g[fAt  depend  upon  the  solution  variables 
For  an  adiabatic  analysis,  z[fAt  and  gf^At  depend  upon  0 ,  via  Eq.  (38)  and  the  temperature  dependence  of  elas¬ 
tic  moduli.  The  temperature  rate  for  a  given  increment  spanning  t  and  t  +  At  is  found,  from  Eq.  (33)  with 
r  —  0,  explicitly  in  terms  of  quantities  at  time  t : 


dm 


4  V  +  ?So(se  :  ee)  +  4div(/tVx0)  |  , 
Pc  c  Pc  , 


where 


p= 


p(M  -  e(de^))t 


(71) 


(72) 


meaning  that  1  -  /?  is  the  ratio  of  time  rates  of  stored  energy  to  plastic  work.  The  temperature  at  time  t  +  Ans 
updated  as 


dt+At  —  dt  +  dAt , 


(73) 


and  the  thermal  deformation  at  the  end  of  the  step  is  determined  from 

ff+A,  =  exp(aT01Af)f?.  (74) 

For  each  time  increment,  Eqs.  (70)  are  solved  implicitly  using  values  of  d ,  dt+&t,  and  fet+At  found  from  (71)— (74). 
The  thermoelastic  term  in  Eq.  (33)  can  be  rewritten,  to  first  order,  as  [19] 

pddg^  :  ee  =  9f-lde( se  :  ee)  -  p0/3e(/'fl“1)( se  :  ee) 

*  drl{[MrpgxS + MdV4 +rsgh)mxS + 3 


(75) 
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For  the  pure  W  grains,  possible  cleavage  fractures  are  captured  by  the  model  features  discussed  in  Section 
3.2.  In  the  numerical  implementation,  prior  to  the  initiation  of  damage  in  a  particular  grain,  criterion  Eq.  (42) 
is  checked  at  the  end  of  the  iteration  cycle  (i.e.,  at  time  t  +  At)  in  that  grain,  for  each  of  the  nine  potential 
damage  planes  ({100}  and  {1 10}),  with  initiation  strengths  computed  via  Eq.  (43)  using  updated  local  tem¬ 
perature  values  6t+At.  Upon  initiation  in  a  particular  crystal,  the  rate  of  crack  extension  for  the  subsequent 

is  computed  from  Eq.  (46),  using  values  of  the  applied  deformation  rate  imposed 


time  increment,  h 


t+At 


(constantly)  over  the  time  increment  and  the  projected  nominal  stress  at  the  end  of  the  increment.  In  the  sub¬ 
sequent  cycle,  the  crack  radius  is  updated  simply  as 


r[k) 

rw 


t+At 


-+) 


(A*) 


I 

rw 


t+At 


(76) 


the  instantaneous  crack  opening  magnitude  is  computed  from  the  resolved  normal  stress  from  Eq.  (45),  and 
the  orientation  of  the  crack  opening  vector  in  the  spatial  configuration  is  updated  by  the  elastic  rotation  as 
specified  in  Eq.  (47). 

Upon  computation  of  the  response  (i.e.,  stress  state,  temperature,  and  damage  variables)  of  each  crystal  of 
each  phase,  the  average  nominal  stress  St+At,  temperature  9t+At,  and  intragranular  damage  deformation  F^| 
are  calculated  from  respective  Eqs.  (61),  (19)  and  (44).  Then,  possible  intergranular  fractures  are  addressed  by 
an  implementation  of  the  model  featured  in  Section  3.3.  Recall  that  in  the  present  scheme,  intergranular  frac¬ 
ture  is  restricted  to  a  single  plane  whose  reference  orientation  M  is  determined  from  a  maximum  principal 
stress-based  criterion.  Before  damage  develops,  the  orientation  of  the  potential  plane,  M,  is  found  as  the  prin¬ 
cipal  eigenvector  of  the  symmetric  nominal  stress  tensor  SAaF~lB  \t+At-  Criterion  Eq.  (50)  is  then  checked  versus 
the  resolved  nominal  stress  due  to  S,+A*,  with  the  initiation  strength  computed  in  Eq.  (51)  using  the  updated 
temperature  value  9t+At •  Upon  initiation  of  intergranular  damage,  the  rate  of  crack  extension  for  the  subse¬ 


quent  time  increment, 


,  is  computed  from  Eq.  (55),  and  in  the  subsequent  cycle,  the  crack  radius  is 

t+At 


updated  in  a  manner  analogous  to  (76).  The  instantaneous  crack  opening  magnitude  is  computed  from  the 
resolved  normal  stress  using  Eq.  (53),  and  the  orientation  of  the  crack  opening  vector  in  the  spatial  configu¬ 
ration  is  updated  by  the  net  material  rotation  as  specified  by  Eq.  (52). 

At  the  conclusion  of  each  iteration  cycle,  the  contributions  from  material  deformation,  W  cleavage,  and 
intergranular  fracture  are  summed  using  Eq.  (62)  to  yield  the  total  deformation  F?+A?.  The  porosity  and  micro¬ 
crack  densities  are  computed  via  Eqs.  (66)  and  (67),  and  then  substituted  into  Eq.  (65)  to  compute  the  damage 
variable  At+At.  Lastly,  the  effective  stress  St+At  is  found  from  Eq.  (63). 

Note  that  the  present  scheme  is  easy  to  implement  if  numerical  crystal  plasticity  routines  are  available,  as 
the  damage  computations  are  effectively  uncoupled  from  the  constitutive  update  of  the  elastoplastic  response 
of  each  crystal.  However,  we  found  that  very  small  integration  time  steps  were  required  to  capture  the  intri¬ 
cacies  of  damage  evolution,  due  to  the  explicit  time  integration  of  the  microcrack  radii,  Eq.  (76).  It  is  under¬ 
stood  that  should  the  model  be  employed  in  large  scale  calculations  (i.e.,  a  setting  with  numerous  finite 
elements,  for  example)  the  plasticity  and  damage  algorithms  may  be  more  efficiently  packaged  in  a  fully  cou¬ 
pled  implicit  scheme,  though  this  will  require  substantial  algorithm  development  due  to  the  unique  nature  of 
the  kinematic  decomposition  Eq.  (62),  and  convergence  difficulties  may  arise  due  to  global  strain  softening  in 
the  damaged  regime.  Furthermore,  as  the  total  deformation  F  cannot  be  an  outcome  in  a  general  large  scale 
simulation,  F  rather  than  F  will  need  to  be  imposed  at  each  integration  point. 


5.  Results  and  discussion 

Predictions  of  the  model  under  uniaxial  tensile  states  of  stress  are  discussed  first.  Deformations  and  strain 
rates  were  prescribed  in  the  simulations  according  to  F 33  =  1  +  st,  with  s  the  nominal  strain  rate  in  the  axial 
direction  (Cartesian  coordinates  are  used  here).  In  the  simulations,  isothermal  conditions  ( 6  =  300  K)  were 
prescribed  for  low  strain  rate  tests  (e  =  0.1/s  and  s  =  0.0001  / s),  and  adiabatic  conditions  were  invoked  at  a 
higher  rate  (e  =  750/s),  with  an  initial  temperature  set  to  90  =  300  K.  As  reported  in  [25],  the  experimental 
data  were  obtained  for  a  93W-5Ni-2Fe  tungsten  alloy  at  low  nominal  strain  rates  (e  =  0. 1/s  and 
e  =  0.0001/s)  with  an  Instron  test  machine,  and  at  the  high  strain  rate  (s  =  750/s)  with  a  split  Hopkinson 
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Fig.  5.  Tensile  stress  versus  tensile  strain. 


bar  apparatus  (i.e.,  Kolsky  bar).  The  experimental  data  are  representative  of  behavior  observed  over  numer¬ 
ous  tests  at  each  strain  rate. 

Numerical  results  from  the  model  are  compared  with  experimental  stress-strain  data  in  Fig.  5.  In  these  sim¬ 
ulations,  random  initial  lattice  orientations  were  assigned  to  the  polycrystalline  aggregate,  consisting  of  300  W 
crystals  and  100  binder  crystals.  In  the  model  and  in  the  experiments,  peak  stresses  (i.e.,  ultimate  tensile 
strengths)  increase  with  increasing  applied  strain  rate,  in  the  former  captured  by  the  strain  rate  sensitivity 
entering  the  flow  rule  Eq.  (37).  Ductility  (i.e.,  elongation  at  failure  or  rupture)  tends  to  decrease  with  increas¬ 
ing  strain  rate,  as  higher  tensile  stresses  more  quickly  activate  a  greater  number  of  damage  sites  (i.e.,  micro¬ 
fractures)  in  the  material.  At  the  lowest  strain  rate,  s  =  0.0001  / s,  the  peak  average  stress  S33  in  the  model, 
1.12  GPa,  is  insufficient  to  activate  intergranular  damage  mechanisms  (strength  of  1.20  GPa  in  criterion 
(50))  or  cleavage  fracture  of  any  W  crystals  (strengths  of  1.75  GPa  and  1.90  GPa  for  {100}  and  {110}  planes, 
respectively).  At  the  intermediate  strain  rate  s  =  0.1/s,  intergranular  mechanisms  are  activated  but  cleavage 
fractures  are  not.  This  prediction  (and  the  assigned  material  properties  in  Tables  2  and  3)  agrees  with  exper¬ 
imental  observations  that  intergranular  failure  modes  (most  notably  W-W  boundaries  which  serve  as  the 
weakest  link  in  the  micro  structure)  tend  to  initiate  more  readily  than  cleavage  fractures,  especially  at  room 
temperature  [24,25,29].  Finally,  at  the  highest  strain  rate  considered,  s  =  750/s,  both  inter-  and  intragranular 
damage  modes  are  active  due  to  the  relatively  high  tensile  stresses,  with  peak  values  of  S 33  on  the  order  of 
1.65  GPa  achieved  in  both  simulation  and  experiment.  At  the  highest  strain  rate,  total  failure  of  the  poly  crys¬ 
talline  aggregate  was  achieved  in  the  simulation  when  the  density  of  cleavage  microcracks  reached  the  limit 
prescribed  in  Table  4,  co^7^ax  =  0.38,  such  that  the  nominal  extension  at  rupture,  F 33  =  1.06,  matched 
the  experimental  result.  Notice  also  in  Table  4  that  critical  values  for  total  rupture  via  intergranular  failure 
/ co™x)  and  porosity  (<pfai1)  are  not  assigned,  as  aforementioned  experiments  as  well  as  numerical  investi¬ 
gations  at  the  microscopic  scale  [19]  have  indicated  that  while  evolution  of  damage  typically  first  occurs  at 
intergranular  sites,  transgranular  fractures  must  be  induced  in  order  to  propagate  a  microcrack  across  a  poly¬ 
crystalline  volume  of  significant  dimensions,  to  a  size  large  enough  to  cause  catastrophic  rupture  of  the  spec¬ 
imen,  as  not  enough  intergranular  sites  exist  in  the  material  to  exclusively  support  such  a  large  crack. 

Figs.  6(a)  and  (b)  depict  the  model  predictions  of  crack  density  and  porosity,  respectively,  at  strain  rates  of 
e  =  0.1/s  and  8  =  750/s.  Damage  initiates  early  in  each  case,  at  an  applied  stretch  around  F 33  =  1.003,  in  rea¬ 
sonable  agreement  with  numerical  results  of  a  previous  micromechanical  investigation  [19].  From  Fig.  6(a),  at 
the  lower  strain  rate  of  £  =  0.1/s,  the  intergranular  crack  density  (Eq.  (67))  saturates  to  its  maximum  allowed 
value  at  F 33  =  1.05,  and  cleavage  fractures  do  not  occur.  On  the  other  hand,  at  the  higher  strain  rate  of 
s  =  750/s,  both  grain  boundary  and  cleavage  fractures  increase  steadily  until  final  rupture  occurs  at 
F 33  =  1.06.  Similar  trends  arise  for  the  porosity  (p  in  Fig.  6(b):  saturation  of  the  volume  fraction  of  intergran¬ 
ular  damage  occurs  at  the  lower  strain  rate  e  =  0.1/s,  whereas  both  inter-  and  intragranular  damage  contrib¬ 
ute  steadily  to  cp  at  the  higher  rate  8  =  750/s  until  rupture.  Notice  that  at  the  lower  strain  rate,  tp  decreases 
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(b) 


Fig.  6.  (a)  Normalized  damage  versus  tensile  strain  and  (b)  damage  versus  tensile  strain. 


very  slightly  with  increasing  deformation  upon  reaching  its  peak  value  (on  the  order  of  (p  =  0.49),  as  some 
elastic  microcrack  closure  takes  place  due  to  stress  relaxation. 

Model  results  for  60  randomly  oriented  crystals  (50  W  grains  and  10  binder  grains)  are  compared  with 
results  of  400  randomly  oriented  crystals  (300  W  grains  and  100  binder  grains)  in  Fig.  7.  It  was  found  that 
including  more  than  400  crystals  resulted  in  negligible  changes  to  the  homogenized  model  predictions  for 
stress-deformation-failure  behavior.  However,  as  is  clear  from  Fig.  7,  differences  in  results  are  apparent  when 
as  few  as  60  grains  are  used.  In  large  scale  computations,  one  may  choose  to  include  fewer  grains  (e.g.,  60  as 
opposed  to  400)  to  reduce  execution  times,  should  discrepancies  on  the  order  of  those  in  Fig.  7  be  deemed 
acceptable. 

Effects  of  various  volume  fractions  of  W  phase  are  shown  in  Fig.  8.  As  the  volume  fraction  of  the  stiff  W 
phase  increases  relative  to  that  of  the  more  ductile  binder,  the  stress  S 33  supported  by  the  aggregate  increases 
at  low  and  high  strain  rates,  in  general  agreement  with  experimental  trends  reported  in  [26].  However,  exper¬ 
imental  data  [26]  also  indicate  that  ductility  (i.e.,  strain  at  dynamic  rupture)  should  decrease  as  the  volume 
fraction  of  W  is  increased,  a  phenomenon  not  captured  by  the  model  in  Fig.  8,  where  at  the  high  strain  rate, 
s  =  750/s,  failure  occurs  consistently  at  F 33  =  1.06  regardless  of  volume  fraction.  A  reduction  in  ductility  with 
increasing  volume  fraction  could  easily  and  logically  be  captured  by  allowing  critical  microcrack  densities  and 
porosity  listed  in  Table  4  to  depend  upon  the  volume  fraction  of  each  phase,  though  more  experimental  data  is 
needed  to  justify  for  these  parameters  particular  choices  of  functions  of  the  initial  micro  structure. 

Finally,  Fig.  9  compares  model  predictions  with  experimental  results  under  torsional  (i.e.,  pure  shear) 
deformation.  Deformations  and  strain  rates  were  prescribed  in  the  simulations  according  to  Fn  =  yt,  with 
y  the  nominal  shear  strain  rate.  Isothermal  conditions  (9  =  300  K)  were  prescribed  for  the  low  strain  rate  test 
(y  =  0.0001/s),  and  adiabatic  conditions  were  invoked  at  the  higher  rate  (y  =  600/s),  with  an  initial 
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temperature  of  60  =  300  K.  Again,  400  randomly  oriented  grains  (300  W  and  100  binder)  were  simulated.  The 
experimental  data  were  obtained  for  a  93W-5Ni-2Fe  tungsten  alloy  with  a  split  Hopkinson  bar  apparatus,  as 
discussed  in  [23].  The  macroscopic  shear  stress  l}2  is  reported  in  Fig.  9,  where  Iab  =  J~lFaASAb  is  the  symmet¬ 
ric  Cauchy  stress.  At  the  lower  strain  rate,  y  =  0.0001  / s,  agreement  is  excellent  between  model  and  experi¬ 
ment,  where  significant  damage  does  not  ensue  in  either  case.  At  the  higher  strain  rate,  y  =  600/s,  the 
strength  predicted  by  the  model  is  excessive  relative  to  that  of  the  experiment  (peak  stress  of 
l}2  =  1.30  GPa  versus  0.94  GPa).  While  mode  I  damage  is  captured  by  the  present  implementation  of  the 
model  (and  these  mechanisms  are  activated  in  the  simulation  at  y  =  600/s),  the  model  does  not  capture  physics 
of  failure  reported  in  the  high  rate  torsion  experiments  [23]:  adiabatic  shear  localization  followed  by  tangential 
fracture  (modes  II  and  III). 

6.  Conclusions 

Presented,  upon  invocation  of  homogenization  methods  predicated  upon  the  generalized  Gauss’  theorem,  is 
a  multiscale  framework  for  describing  finite  deformation  and  failure  mechanisms  in  multi-phase  polycrystals. 
Damage  contributions  arising  from  displacement  discontinuities  such  as  cracks,  voids,  and  shear  bands  are 
captured  explicitly  in  the  framework’s  kinematics  and  balance  relations  through  additive  decompositions  of 
the  total  deformation  gradient  (Eqs.  (5)  and  (62))  and  nominal  stress  (Eqs.  (11)  and  (61)).  The  present 
approach  allows  for  both  strain  and  rotation  accommodation  due  to  damage  evolution. 

A  specific  version  of  the  framework  has  been  implemented  to  describe  the  thermomechanical  response  of  a 
two-phase  tungsten  alloy.  The  complete  implementation  consists  of  single  crystal  plasticity  models  capturing 
thermo-elasto-viscoplasticity  in  each  phase,  a  cleavage  fracture  model  representing  failure  of  tungsten  grains, 
an  intergranular  fracture  model  describing  damage  evolution  at  grain  and  phase  boundaries,  a  polycrystal 
averaging  scheme  with  Taylor  constraints,  and  a  damage  homogenization  model  relating  the  effective  stress 
to  porosity  and  crack  densities.  Orientation-dependent  damage  initiation  is  specified  on  preferred  cleavage 
planes,  with  relative  fracture  strengths  scaled  in  agreement  with  reported  trends  in  fracture  toughness,  and 
ductile-to -brittle  transitions  are  captured  to  first  order  by  temperature-dependent  initiation  criteria.  Model 
predictions  for  stress-versus-deformation  behavior  agree  with  experimental  results  at  low  and  high  tensile 
strain  rates  and  low  shear  strain  rates.  The  model  is  unable  to  match  experimental  high  rate  shear  data,  as 
sub-models  for  shear  localization  and  tangential  cracking  are  not  included  in  the  present  implementation, 
and  remain  to  be  developed  in  future  studies.  Additional  areas  for  model  development  and  improvement, 
as  more  experimental  data  become  available  for  validation,  include  development  of  a  fully  anisotropic  damage 
effect  tensor  entering  Eq.  (63),  relaxation  of  the  Taylor  deformation  gradient  restriction  Eq.  (60)  to  capture 
intergranular  interactions,  numerical  incorporation  of  heat  conduction,  and  further  exploration  of  the  homo¬ 
genized  thermodynamic  framework  developed  in  Eqs.  (19)-(24)  and  Eqs.  (68)  and  (69). 

Acknowledgement 

This  work  was  supported  by  the  US  Army  Research  Laboratory. 

References 

[1]  D.  Krajcinovic,  Damage  Mechanics,  North-Holland,  Elsevier,  Amsterdam,  1996. 

[2]  A.L.  Gurson,  Continuum  theory  of  ductile  rupture  by  void  nucleation  and  growth:  Part  1 — yield  criteria  and  flow  rules  for  porous 
ductile  media,  J.  Eng.  Mater.  Sci.  Technol.  99  (1977)  2-15. 

[3]  D.J.  Bammann,  E.C.  Aifantis,  A  damage  model  for  ductile  metals,  Nucl.  Eng.  Des.  116  (1989)  355-362. 

[4]  G.R.  Johnson,  W.H.  Cook,  Fracture  characteristics  of  three  metals  subjected  to  various  strains,  strain  rates,  temperatures,  and 
pressures,  Eng.  Fract.  Mech.  21  (1985)  31^18. 

[5]  S.  Murakami,  Notion  of  continuum  damage  mechanics  and  its  application  to  anisotropic  creep  damage  theory,  J.  Eng.  Mater. 
Technol.  105  (1983)  99-105. 

[6]  A.  Menzel,  P.  Steinmann,  Geometrically  non-linear  anisotropic  inelasticity  based  on  fictitious  configurations:  application  to  the 
coupling  of  continuum  damage  and  multiplicative  elasto-plasticity,  Int.  J.  Numer.  Methods  Eng.  56  (2003)  2233-2266. 

[7]  R.S.  Kumar,  R.  Talreja,  A  continuum  damage  model  for  linear  viscoelastic  materials,  Mech.  Mater.  35  (2003)  463-480. 


J.D.  Clayton  /  Theoretical  and  Applied  Fracture  Mechanics  45  (2006)  163-185 


185 


[8]  A.M.  Rajendran,  D.J.  Grove,  Modeling  the  shock  response  of  silicon  carbide,  boron  carbide  and  titanium  diboride,  Int.  J.  Impact 
Eng.  18  (1996)  611-631. 

[9]  H.D.  Espinosa,  P.D.  Zavattieri,  S.K.  Dwivedi,  A  finite  deformation  continuum  discrete  model  for  the  description  of  fragmentation 
and  damage  in  brittle  materials,  J.  Mech.  Phys.  Solids  46  (1998)  1909-1942. 

[10]  G.Z.  Voyiadjis,  R.K.  Abu  Al-Rub,  A.N.  Palazotto,  Thermodynamic  framework  for  coupling  of  non-local  viscoplasticity  and  non¬ 
local  anisotropic  viscodamage  for  dynamic  localization  problems  using  gradient  theory,  Int.  J.  Plasticity  20  (2004)  981-1038. 

[11]  T.E.  Lacy,  D.L.  McDowell,  R.  Talreja,  Gradient  concepts  for  evolution  of  damage,  Mech.  Mater.  31  (1999)  831-860. 

[12]  S.  Ahzi,  S.  Schoenfeld,  Mechanics  of  porous  polycrystals:  a  fully  anisotropic  flow  potential,  Int.  J.  Plasticity  14  (1998)  829-839. 

[13]  A.  Zubelewicz,  Micromechanical  study  of  ductile  polycrystalline  metals,  J.  Mech.  Phys.  Solids  41  (1993)  1711-1722. 

[14]  H.W.  Meyer,  D.S.  Kleponis,  Modeling  the  high  strain  rate  behavior  of  titanium  undergoing  ballistic  impact  and  penetration,  Int.  J. 
Impact  Eng.  26  (2001)  509-521. 

[15]  T.  Watanabe,  Grain  boundary  design  for  the  control  of  intergranular  fracture,  Mater.  Sci.  Forum  46  (1989)  25-48. 

[16]  J.D.  Clayton,  D.L.  McDowell,  Finite  polycrystalline  elastoplasticity  and  damage:  multiscale  kinematics,  Int.  J.  Solids  Struct.  40 
(2003)  5669-5688. 

[17]  C.  Teodosiu,  F.  Sidoroff,  A  finite  theory  of  elastoplasticity  of  single  crystals,  Int.  J.  Eng.  Sci.  14  (1976)  713-723. 

[18]  P.  Gumbsch,  J.  Riedle,  A.  Hartmaier,  H.F.  Fischmeister,  Controlling  factors  for  the  brittle-to-ductile  transition  in  tungsten  single 
crystals,  Science  282  (1998)  1293-1295. 

[19]  J.D.  Clayton,  Dynamic  plasticity  and  fracture  in  high  density  polycrystals:  constitutive  modeling  and  numerical  simulation,  J.  Mech. 
Phys.  Solids  53  (2005)  261-301. 

[20]  R.  Hill,  On  constitutive  macro-variables  for  heterogeneous  solids  at  finite  strain,  Proc.  R.  Soc.  Lond.  A326  (1972)  131-147. 

[21]  S.  Nemat-Nasser,  Averaging  theorems  in  finite  deformation  plasticity,  Mech.  Mater.  31  (1999)  493-523. 

[22]  G.I.  Taylor,  Plastic  strain  in  metals,  J.  Inst.  Metals  62  (1938)  307. 

[23]  T.  Weerasooriya,  P.A.  Beaulieu,  Effects  of  strain  rate  on  the  deformation  and  failure  behavior  of  93W-5Ni-2Fe  under  shear  loading, 
Mater.  Sci.  Eng.  A  172  (1993)  71-78. 

[24]  T.  Weerasooriya,  P.  Moy,  R.  Dowding,  Effect  of  W-W  contiguity  on  the  high  shear  strain  rate  behavior  of  93W-5Ni-2Fe  tungsten 
heavy  alloy,  in:  A.  Bose,  R.  Dowding  (Eds.),  Proceedings  of  the  2nd  International  Conference  on  Tungsten  and  Refractory  Metals, 
Metal  Powder  Industries,  1994,  pp.  401-409. 

[25]  T.  Weerasooriya,  Deformation  and  failure  behavior  of  a  tungsten  heavy  alloy  under  tensile  loading  at  different  strain  rates, 
in:  Proceedings  of  the  SEM  Annual  Conference  on  Experimental  Mechanics,  Charlotte,  NC,  USA,  2-4  June  2003. 

[26]  T.  Weerasooriya,  P.  Moy,  High  shear  strain  rate  behavior  of  W-Ni-Fe  tungsten  heavy  alloy  composites  as  a  function  of  matrix 
volume  fraction,  ARL-TR-1964,  1998. 

[27]  S.E.  Schoenfeld,  D.J.  Benson,  Modeling  penetration  mechanisms  and  ballistic  performance  in  high  aspect  ratio  tungsten  single-crystal 
rods:  a  crystal  plasticity  model  suitable  for  impact  calculations,  in:  S.N.  Atluri,  G.  Yagawa  (Eds.),  Advances  in  Computational 
Engineering  Science,  Tech  Science  Press,  1997,  pp.  1116-1121. 

[28]  Z.  Baoping,  Y.  Chao,  X.  Yingming,  J.  Chunlan,  Responsive  behaviour  of  93  wt.%  tungsten  alloy  under  the  intense  shock  loading,  in: 
Z.  Zhemin,  T.  Quingming  (Eds.),  Proceedings  of  IUTAM  Symposium  on  Impact  Dynamics,  Peking  University  Press,  Peking,  PRC, 
1994,  pp.  283-297. 

[29]  R.G.  O’Donnell,  R.L.  Woodward,  Influence  of  temperature  on  the  fracture  of  a  W-Ni-Fe  alloy,  J.  Mater.  Sci.  35  (2000)  4319-4324. 

[30]  A.  Molinari,  G.R.  Canova,  S.  Ahzi,  A  self-consistent  approach  of  the  large  deformation  polycrystal  plasticity,  Acta  Metall.  35  (1987) 
2983-2994. 


NO.  OF 

COPIES  ORGANIZATION 


1  DEFENSE  TECHNICAL 
(PDF  INFORMATION  CTR 
ONLY)  DTICOCA 

8725  JOHN  J  KINGMAN  RD 
STE  0944 

FORT  BEL  VOIR  VA  22060-6218 

1  US  ARMY  RSRCH  DEV  & 
ENGRG  CMD 
SYSTEMS  OF  SYSTEMS 
INTEGRATION 
AMSRD  SS  T 
6000  6TH  ST  STE  100 
FORT  BEL  VOIR  VA  22060-5608 

1  DIRECTOR 

US  ARMY  RESEARCH  LAB 
IMNE  ALC  IMS 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1 197 

3  DIRECTOR 

US  ARMY  RESEARCH  LAB 
AMSRD  ARL  Cl  OK  TL 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1 197 


ABERDEEN  PROVING  GROUND 
1  DIR  USARL 

AMSRD  ARL  Cl  OK  TP  (BLDG  4600) 


NO.  OF 

COPIES  ORGANIZATION 


ABERDEEN  PROVING  GROUND 


32  DIR  USARL 

AMSRD  ARL  Cl  HC 
P CHUNG 
AMSRD  ARL  WM 
J  MCCAULEY 
T  WRIGHT 

AMSRD  ARL  WM  MA 
W  NOTHWANG 
AMSRD  ARL  WM  TA 
S  SCHOENFELD 
AMSRD  ARL  WM  TC 
M  FERMEN  COKER 
R  COATES 
AMSRD  ARL  WM  TD 
S  BILYK 
T  BJERKE 
D  CASEM 
J  CLAYTON  (5  CPS) 
T  CLINE 
D  DANDEKAR 
W  EDMANSON 
M  GREENFIELD 
C  GUNNARSSON 
Y  HUANG 
K  IYER 
R  KRAFT 
BLOVE 
S  MCNEILL 
H  MEYER 
RMUDD 
M  RAFTENBERG 
E  RAPACKI 
M  SCHEIDLER 
S  SEGLETES 
T  WEERASOORIYA 


Intentionally  left  blank. 


